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Messenger RNA translation is often studied by means of statistical-mechanical models based on 
the Asymmetric Simple Exclusion Process (ASEP), which considers hopping particles (the ribo- 
somes) on a lattice (the polynucleotide chain). In this work we extend this class of models and 
consider the two fundamental steps of the ribosome's biochemical cycle following a coarse-grained 
perspective. In order to achieve a better understanding of the underlying biological processes and 
compare the theoretical predictions with experimental results, we provide a description lying be- 
tween the minimal ASEP-like models and the more detailed models, which are analytically hard 
to treat. We use a mean-field approach to study the dynamics of particles associated with an in- 
ternal stepping cycle. In this framework it is possible to characterize analytically different phases 
of the system (high density, low density or maximal current phase). Crucially, we show that the 
transitions between these different phases occur at different parameter values than the equivalent 
transitions in a standard ASEP, indicating the importance of including the two fundamental steps 
of the ribosome's biochemical cycle into the model. 

PACS numbers: 87.10.-e, 05.70.Ln, 87.16.A- 



I. INTRODUCTION 

The translation of the messenger RNA (mRNA) is the 
final step of protein synthesis. During this process the 
information enclosed in the triplet code of the nucleotide 
chain is translated into the amino acid sequence of the en- 
coded proteins. Translation is usually viewed as a three- 
stage process [T]-[3]: during initiation a ribosome (com- 
plex of proteins and RNA) binds the mRNA molecule 
(a sequence of nucleotides previously transcribed from 
the DNA) and after a series of biochemical reactions, it 
moves along the chain. This stage in which the protein 
is built up amino acid by amino acid according to the 
mRNA sequence is called elongation. Each elongation 
step consists in turn of a series of biochemical reactions 
which define the ribosome's biochemical cycle. Lastly, 
the ribosome reaches the termination codon and leaves 
the mRNA releasing the protein. This last step is called 
termination. In this paper we propose a model for the 
elongation stage of mRNA translation. 

This process, primarily controlled by the dynamics of 
ribosomes along the mRNA chain, bears a resemblance 
to a one-dimensional driven lattice gas. For this reason, 
the mRNA translation inspired a statistical-mechanical 
class of models known as Asymmetric Simple Exclusion 
Processes (ASEPs). They have been introduced in the 
biophysical literature as models representing the dynam- 
ics of ribosomes along an mRNA chain [4j [5] . Later, this 
class of model has been studied from a more theoretical 
point of view [rjHTB] and the possible biological applica- 
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tions have been rediscovered only recently, not only for 
protein synthesis [12 [TV 20J but also for the movement 
of molecular motors [2P - f23] . Other non-biological appli- 
cations have been studied too (see, e.g., [H]). A detailed 
discussion of the ASEP can be found, for instance, in 

Refs. num. 

A traditional exclusion process consists of particles 
moving along a lattice with only steric interactions. In 
other words, each lattice site can be occupied just by one 
particle at a time. Although this approach is interest- 
ing from a theoretical point of view, it is however not 
a realistic way of describing translation since it encom- 
passes the whole ribosomal elongation cycle in a single 
step. Other approaches to modelling translation include 
several (up to fifteen) distinct phases of the ribosome's 
mechano-chemical cycle [251427J . We want to place our- 
selves between the class of minimal models and the more 
detailed models, which are difficult to analyze. Conse- 
quently, we consider two fundamental steps of the ribo- 
some's biochemical cycle following a coarse-grained pic- 
ture (Section III I. We show that former ASEP-like mod- 



els correspond to a limiting case of the model we use. As 
will become clear from the discussion, this limiting case 
is however biologically not plausible, showing the need 
for an extension of previous models. 

In the following section we explain the biological frame- 
work (in particular the role of ribosomes and transfer 
RNAs) and then in Section III we introduce the model 
from a mathematical point of view. The results for pe- 
riodic and open-boundary systems are presented in Sec- 
tion IV Using the same approach as Ref. [5] we analyze 



the model in the mean-field approximation and compare 
this model with a typical exclusion process. In addition, 
we show that the same results can be achieved by us- 
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ing an extremal principle [25U30] , Although we observe 
the same variety of phase transitions that one would ex- 
pect from the simpler case without the particle's internal 
states, the locations of the critical points change substan- 
tially and depend on the internal dynamics of the parti- 
cles. Finally, in Section [V] we discuss the results and the 
conclusions from a theoretical and biological viewpoint. 

II. BIOLOGICAL BACKGROUND 

Here we briefly introduce the underlying biological pro- 
cess that we want to describe. More information can be 
found, e.g., in Refs. [THS]. 

The mRNA is a nucleotide chain composed of four dif- 
ferent bases (A,C,G,U); each group of three nucleotides 
is called a codon and specifies a certain amino acid. The 
keys for deciphering the code (the sequence of codons) are 
the transfer RNAs (tRNAs), freely diffusing molecules 
carrying amino acids. They have a region composed 
of three nucleotides (the anticodon) matching the cor- 
responding codon on the mRNA. Moreover, tRNAs with 
the same anticodon transport the same amino acid. In 
most species, there are 35-40 distinct species of tRNA, 
each responsible for decoding a particular set of codons. 
Ribosomes are complexes of proteins and RNA which 
move along the mRNA in a fixed direction (from the 5' 
to the 3' region, i.e. from the part of the chain that has 
been transcribed first towards the other end) and cat- 
alyze the assembly of amino acids delivered by tRNAs. 
Each ribosome has three regions of interaction for the 
tRNA. They are called Aminoacyl (A), Peptidyl (P) and 
Exit (E) sites (see Fig. [l}. 




- codon 



FIG. 1: (a) Illustration of a ribosome with empty A, P, and 
E sites along an mRNA chain, (b) Classical representation 
of a tRNA with its anticodon in the lower region and the 
corresponding amino acid (aa) bound together. 

The main steps of the elongation process are shown in 
Fig. [2] Following the translation initiation, a tRNA is in 
the P site bound with the first amino acid of the grow- 
ing polypeptide chain. Then, a complex of EFlev • GTP 
• tRNA aa (with tRNA aa we denote a tRNA bound to 
the amino acid aa) diffuses into the empty A site on the 
codon at position i. If the anticodon of the tRNA cannot 
base-pair with the codon, then the complex EF1 • GTP 



• tRNAaa is released and the process is repeated until a 
correct tRNA binds the ribosome. If on the other hand 
the correct tRNA anticodon base-pairs with the corre- 
sponding codon on the mRNA, then GTP is hydrolyzed 
and a conformational change in the ribosomal structure 
occurs. This change leads to the transfer of the nascent 
peptide from the P-site tRNA, to the amino acid car- 
ried by the A-site tRNA. The altered structure of the 
ribosome does not allow the cognate tRNA to unbind 
and leave the chain. Following the peptidyl transfer re- 
action and the incorporation of the new amino acid into 
the growing polypeptide chain, the ribosome translocates 
one codon (assuming ribosome progress is not blocked 
by any stalled ribosomes at downstream positions on the 
mRNA). The translocation process is catalyzed by the 
complex EF2 • GTP which induces a second conforma- 
tional change in the ribosome. The tRNA at the P site 
is then transferred to the E site. The ribosome is now 
back to the first step of its biochemical cycle, with the 
growing polypeptide chain bound to the tRNA in the P 
site and the empty A site on the codon i + 1, ready to 
receive another tRNA complex. This process is iterated 
until the end of the mRNA chain, where the ribosome dis- 
associates from the system and releases the synthesized 
protein. As soon as initiating ribosomes have moved suf- 
ficiently downstream to create space at the beginning of 
the mRNA, a new ribosome can bind the polynucleotide 
chain. Thus, several ribosomes can translate the same 
mRNA at the same time. 




FIG. 2: Sketch of the translation elongation process. Once 
the ribosome finds the cognate tRNA (1), the tRNA in the E 
site abandons the ribosome (2) and the peptide carried by the 
existing, P-site tRNA binds the amino acid on the new, A-site 
tRNA (3). At that point the ribosome translocates provided 
that the next codon is empty (4). The ribosome is thus in the 
position to accept another tRNA and iterate the elongation 
till the end of the mRNA chain. 



Experimental data strongly indicate that searching 
for the correct tRNA, and not the translocation, is the 
rate limiting step of the biochemical cycle of the ribo- 
some [31] [32]. Therefore, in this work we shall approx- 
imate the whole biochemical cycle of the ribosome by a 
two-state cycle: (i) searching for the correct tRNA and 
(ii) translocation from one codon to the next. 
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From the modelling point of view, we shall consider 
particles changing an internal degree of freedom, or state, 
which influences their motion. Thus, a ribosome is rep- 
resented by a two-state particle denoting the absence of 
the cognate tRNA in its Aminoacyl site (state 1) or its 
presence (state 2) as outlined in Fig. [3] 

Throughout this work we assume that the concentra- 
tion of charged tRNAs is homogeneous and large enough 
to neglect fluctuations (the transition rates do not change 
with time). The effects of limited resources [33] and their 
3D diffusion in the cytoplasm [M] are not taken into con- 
sideration here. The individual charged tRNA concentra- 
tions govern the transition rates, which in general depend 
on the particular type of codon. 



correct tRNA with 
rate k 



translocation with 
rate y 




FIG. 3: The two states that a particle assumes represent a 
ribosome waiting for the cognate tRNA (state 1) and ready 
to translocate (state 2). The transition from "state 1" to 
"state 2" occurs with rate k (in general depending on the 
codon i) which mainly models the concentration of tRNA aa . 
The translocation occurs with rate 7. The transitions are not 
reversible. 



III. THE MODEL 

We describe the mRNA molecule as a lattice of discrete 
sites, each one representing one codon. Ribosomes are 
represented by particles hopping from one site of the lat- 
tice to the next. From a more mathematical perspective, 
the occupation number = 0, 1, 2 of the site i describes 
the different states in which it can be found. We say that 
a site is empty if its occupation number is 0. A given site 
i occupied by a particle in state 1 or 2 is respectively 
described by rij = 1 or rij = 2. The set 77 = {ni, . . . , nj,}, 
where L is the length of the lattice, will give the config- 
uration of the system. The only transitions allowed are 
the following: 



1 ->■ 2 
20 ->■ 01 



with rate fc, 
with rate 7, 



(la) 
(lb) 



where the first line means that a particle in state 1 at 
site i changes into state 2 with rate fej, which in general 
depends on the site i. The second line is a schematic 
representation of the translocation of a particle in state 2 



to the next site. Notice that a hopping particle is carried 
back to the state 1. The values of n^'s change with time 
according to these dynamical rules. This dynamics has 
been first introduced in the literature by Klumpp and 
coworkers in [35] to model the traffic of molecular motors 
on a filament. 

The mean density pi of particles on site i can be written 
discerning the contribution of particles in state 1 and 
particles in state 2. Thus, we shall say that A^ is the 
mean density of particles in state 1 at site i and <ii is the 
analogous for particles in state 2. One can write these 
densities in terms of the occupation numbers: 



(n»(2 - n,-)) 
Tk(ni - 1) 



A; 



The mean density of particles at site i is then given by 
Pi = Xi + Gi. The brackets indicate the average of the 
quantities over time. 

The lattice in consideration may have periodic or open 
boundary conditions. We shall study these cases in Sec- 
tions |IV A| and |IVB[ respectively. 

From now on we study the case in which fc, = k Vi 
and use the mean-field approximation, i.e. we neglect 
correlations between the sites ((r^n^) ~ (ni){rij)). These 
approximations simplify the analysis considerably and, as 
we shall show later, yield qualitatively the same results. 

With these prescriptions, the mean-field equations de- 
scribing the evolution of the densities at site % read as 
follows: 



~dt 
dcr l 
~dt 



<7i_i(l - Xi - Oi)7 
k\i — aAl — A i+ i - 



- k\i 



(2a) 
(2b) 



The current ,P + ( ,P_) is defined as the number of particles 
entering (leaving) the site i per unit time. We can write 
the expression of the incoming and outgoing currents by 
using Eqs. §: 



J\ = a 

ji 



i(l - A,; - «7j)7 
cr,(l - Aj+i - <r l+ i)-i 



In this work we consider the steady-state condition 
( <ft = * = Vi), where the currents are the same 
along the lattice (J := J l _ = J'^ Vi). 



IV. RESULTS 

First we study the effects of two-state particles in 
closed lattices and then we investigate the boundary- 
induced phase transitions in open systems. The pre- 
dictions of the mean-field theory are then compared to 
numerical simulations performed with a Bortz-Kalos- 
Lebowitz-like algorithm |36j modified for the dynamic 
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rules ([I]), i.e. a continuous-time Monte Carlo which uses 
a random sequential updating scheme. The first 10 6 iter- 
ations of the algorithm are disregarded. Then, with the 
system in the steady-state, data is collected every 100 
iterations, for a total number of 10 6 iterations. 



A. Periodic-boundary conditions 

Since all sites in a lattice with periodic-boundary con- 
ditions are identical (in the special case of fc$ = k Vz), we 
write Eqs. ^ without the indices i. Note that the mean 
number of particles in state 1 is equal to the local den- 

= A. The same holds for 

L 



sity A, i.e., L 1 {Y, i=1 M n , 

particles in the upper state: (2L)~ 1 (Y^^Li ni$m,2 
One obtains the following equations for the current of 
particles J and for the densities: 



a . 



J = a(l - A - cr)7 



A = 



J 

k 

p-X 



j 
k 



(3) 

(4a) 
(4b) 



The density p plays the role of the control parameter. 
For this reason, it is useful to write Eqs. ^ and Q as 
follows 



J = 



pQ- - p) 



+ (i-p) 



(5) 



(6) 



These results have been previously obtained in [35j . Note 
that both densities A and a are functions of the ratio k/j. 
This result is not an artifact of the mean-field approxi- 
mation, since simulations confirm this dependence (data 
not shown). For the sake of simplicity, in simulations we 
can therefore set the value of 7 to a fixed value (for in- 
stance 7 = 1). From Eq. |5]) one can obtain the value of 
p for which the current is maximal 



1 



k 

7 



1 



2 

k 



i-x. 



(7) 



where \ := (fc/7)(yl + 7/fc — 1), and the maximal value 
of the density of particles in the state 1 (proportional to 
the current J) 



A* := 1 



2k 



(8) 
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FIG. 4: (Color online). Simulations (dashed lines) and mean- 
field (MF) approximation (full lines) of the densities A in black 
and g in green (light gray). This figure shows the curves for 
a ring with N = 250, 7 = 1 and k = 0.05 (a), k = 0.5 (b), 
k = 1 (c), k = 1.5 (d). 



Therefore, 



a := p - A = x 



(9) 



Figure [4] shows the mean-field solutions for the densities 
A and a depending on p, together with the Monte Carlo 
simulations. The overall agreement between the mean- 
field approximation and the simulations is very good. 
Only when fc/7 is very small the discrepancy between the 
mean- field and the simulations becomes large (Fig. [4^,) , 
as already observed in [33]. Thus, in the case fc/7 <C 1, 
correlations are no longer negligible and the mean-field 
overestimates the current J (or, equivalently, the amount 
of particles in the inactive state = 1). However, ana- 
lytical calculations and simulations show the same quali- 
tative behavior and therefore the mean-field approxima- 
tion is sufficient to capture the main features of the un- 
derlying system. 

One can see from Eq. Q and Fig. [5] that as fc/7 in- 
creases, p* approaches the value 0.5 and the current pro- 
file becomes symmetric with respect to the density p. 
This limiting case corresponds to neglecting the internal 
state of the particles, i.e., considering that the transition 
rii = 1 — » rii = 2 occurs instantaneously. Note that in 
this case we recover the results of an ASEP model with 
a single hopping rate 7. 

With increasing p there exist different regimes charac- 
terized by different amounts of particles in state 1 or 2. If 
k < 7, then the curves of A and a cross at pa '■= 1 — k/"f, 
defining two distinct regimes: for p < p^, there is a 
regime in which A > er, i.e., the density of sites with 
rii = 1 is larger than the density of sites with rii — 2. In 
other words, the mRNA is mainly populated by empty 
(vacant A site) ribosomes. For p > p^ we have the op- 
posite situation with A < a, i.e. the ribosomes have the 




FIG. 5: Plot of p* as a function of k/j, Eq. Q. 



tRNA in their A site and are waiting to hop. These dif- 
ferent regimes exist only when k < 7, otherwise A < a 
always. Notice that in general p* ^ pd and therefore, 
this transition is different from the queueing transition. 
These basic observations, though being a simple study of 
the ratio between densities, might reveal an interesting 
biological mechanism (see Section M). 



B. Open-boundary conditions 

In this section we focus on the open-boundary con- 
ditions and present the corresponding results. We first 
discuss the outcomes of the model using an iterative map 
obtained from the mean- field Eqs. §2$. Then we recover 
the same results using an extremal principle. 

A new particle enters the unidimensional lattice with 
rate a representing the translation inititation. After the 
injection of a new particle, the first site is set to have 
m = 1. As usual, the presence/absence of a particle and 
its state are represented by the occupation number Hi 
and in the bulk the dynamics follows the above rules (JlJ . 
Finally, particles abandon the end of the lattice (when 
«i = 2) with probability per unit time (3 (translation 
termination). With these prescriptions it is clear that 
Eqs. ([2]) hold in the bulk, but have to be modified at the 
left boundary (injection) 



~dt 
da 1 
~dt 



a(l - Ai - ax) - kXi (10a) 

= k\ x -ai(l-X 2 -a 2 )7, (10b) 
and at the right boundary (depletion) 

= a L -i(l - X L - a L )^ - kX L (11a) 

kX L - pa L . (lib) 



dXL 
dt 
da L 
~df 



Equations ^ together with the steady-state condition 
lead to the following recursive map for the densities ai 



a i+ i = 1 - J 



1 



1 



k 7<7, 



(12) 




(c) 



FIG. 6: Graphical representation of the recursive map ( 12 1 . 
(a) When J < fc(l — 2\) there are two different fixed points 
(7— , 0+ that collapse when J = k(l — 2\) (b). Panel (c) shows 
the map for the finite size case. 



The fixed points of this map are as follows 
J\ . [7. ~J~' 



<7± 




1 - 



4J 

7 



one of which is stable (0+) and the other unstable (c_). 
In an iterative map like Eq. (12), a^\ is said to be the 



homographic function of ai and, knowing the value of the 
starting point o^, it is possible to find the general term 
ai of the recursion: 



a-a + (a\ 



-a_a+{a l + 2 - <?- 2 ) + a 1 {a i ^ 1 - at 1 ) 



(13) 



Following the approach presented in [5] by Derrida and 
coworkers, we reconstruct the phase diagram of the sys- 
tem by varying the injection rate a and the depletion 
rate j3 (which are both considered to be smaller than 7). 
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FIG. 7: (Color online). Profiles of the density <r,: calculated 
from Eq. (12 1, (a) When the fixed points exist, the density 
can either start close to <r_ and then go away from this value 
at the end of the chain (black -lower- dots) or reach the value 
of the stable fixed point a+ after few iterations [green (light 
gray) dots] . (b) Typical density profile for the MC region with 
high density close to the left boundary and low density close 
to the terminating site. 
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Using Eq. ( 13 1 we can calculate o~l as a function of <j\ 
and J: 



<?L = 0l(<7i, J) 



(14) 



Equation (14 1, together with Eqs. (lOi) and ([Tl~|b) in the 
steady-state condition, determines the values of a±, ul 
and J as a function of a and /3. Moreover, note that 
in the mean-field approximation the densities A^ can be 
readily calculated from A^ = J/k. The system shows 
three different regimes that can be characterised by rea- 



soning on the graphical representation of Eq. ( 12 ) (see 
Fig.§. 

Low Density phase (o~\ ~ <r_, < <j+). If we start to 



iterate the map ( 12 1 close to the unstable fixed point <r_ , 
at the beginning of the lattice the values of o~i remain 
close to this value and then move away (black dots in 
Fig(7^). This is the so-called Low Density phase (LD). 
The recursion (12) and Eqs. (10 1 and (11) provide the 
solutions 



a 

7 



ak(j — a) 
f3j(k + a) 



J 



afc(7 — a) 



7 (fc + a) 

These equations are valid as long as 



IX 



(3>a, 



(15) 



since otherwise Eqs. ( 10 1, ( 11 ) and ( 12 ) are not consistent 
with the conditions <j\ = <r_ and ctl < <x+. The critical 
value a c determines the boundary of the LD regime. 

High Density phase (a\ > cr_, o~l — &+)• Similarly, 
the High Density phase (HD) is reached starting from 
a value a\ > er_. Iterating the map, we reach a value 
<7j arbitrarily close to the stable point o+ [green (light 
gray) dots in Fig(7^]. The initial point a\ lies therefore 
in the domain of attraction of o~ + . Following the previous 
procedure, one obtains the solutions 



(71 = 



akj - (3kj + a(3 2 + k/3 2 



cr/CS + ft) 



o-l = 



J 



7 (fc + /3) 



These solutions exist when 



IX 



< a . 



(16) 



It is worth noting that the critical points a c and f3 c de- 
limiting the LD and the HD phases are functions of k 
and 7. 

Maximal Current phase (o~\ 2 _1 (1 — J/k), ol ^ 
2 _1 (1 — J/k)). This regime is reached when the two 
fixed points collapse and the lattice carries the maximal 
current allowed. This phase occurs when 



P>Pc 



(17) 



and we have the solutions 
a x = 1 



J J 

k a 



ah 



J 



J = fc(l - 2 X ) . 

We expect these results to hold in the limit L —> 00 
and the finite-size effects to be similar to the ones of 
standard ASEP 6J. Thus, in a finite-size system, the 
recursion ( |12[ ) would not have any real fixed points and 
the graphical representation of Fig. [6]d would have to 
be modified into Fig. |6j:. The role of the limited length 
L needs further investigation, but this analysis goes 
beyond the scope of this paper. 



Until now we have used Eq. ( 12 ) as the starting point 



to characterise the different regimes of the process. The 
Maximal Current Principle (MCP) is another viable ap- 
proach which leads to the same results. It was first pre- 
sented by Krug [2S] and then later extended [55] GUI ■ Ac- 
cording to this principle, the boundaries are substituted 
by reservoirs of particles and the dynamics between the 
reservoirs and the lattice is assumed to be the same as 
in the bulk. The MCP states that the current J of an 
open-boundary lattice in the MC regime is given by 



J = max </(p), 

pe[pt+i,p ] 



(18) 



where po and pl+i are respectively the densities of the 
reservoirs of particles at the left and the right boundaries. 
J(p) is the expression of the current as a function of the 
density p that in the bulk we can consider to be given 
by Eq. |5]). The densities po and pl+i are chosen to 
realize the injection and depletion parameters a and /3. 
Equation ( |18[ ) is valid for systems in which the current 
profile has only one maximum, and has to be modified 
if J(p) presents minima [53 ED]. There are no general 
prescriptions for choosing the correct densities po and 
Pl+i °f the reservoirs [T7]. Here we propose a way to 
fix po and pL+i and relate them to the injection and 
depletion rate; with these values we recover the results 
obtained above. 

If we imagine having a reservoir of particles or an extra 
site at i = with density of particles po = Xq + <7o, the 
parameter a can be written as a = jP(no = 2), where 
P(n = 2) is the probability of having the occupation 
number of the site i = equal to 2, i.e., having a parti- 
cle ready to hop from the reservoir to the lattice. Since 
P(no = 2) = (7q one may write [41] : 

a = <7o7 ■ 

On the right boundary we can assume that the density of 
particles at the extra-site L + 1 is related to the depletion 
rate /3 as follows 

/3 = (1 - Pl+i)7 • 
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Now, bearing in mind that x = a * an d using Eq. 0, the 
maximal principle yields the location of the critical points 
by equating cr with a* = p* — X* and pl+i with p* . The 
transitions occur at the same values a c and /3 C obtained 
before in Eqs. (15 1 and (16 1. The current and the bulk 



densities arc then given by the following equations: 



' afc(7 — a) 
•y(k + a) 
f3k(j - p) 



for a < P < 7X (LD) 
for p < a < 7X (HD) 



l(k + P) 

U(1-2 X ) fora,/3> 7X (MC) 



p Q for a < p < 7X (LD) 
/' = <i Pl+i for P < a < 7X (HD) 
1 - X for a, /3 > 7X (MC) 
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FIG. 8: (Color online). Phase diagram of the ASEP with two- 
state particles. Different colored lines correspond to different 
critical values obtained with changing k/j. The MC region is 
larger for low values of k/y and approaches the MC region of 
the standard ASEP in the limit fc/7 — > 00. The dashed line 
separating the LD and HD denotes the first order transition 
between these two regions. 



Now we are finally able to construct the rich phase di- 
agram of the system (Fig. [8]). The model shows the same 
variety of phase transitions of "standard" (particles with- 
out internal states) ASEPs, but the borders between the 
different phases crucially depend on both k and 7. There 
are three different regimes (LD, HD, MC) and the tran- 
sitions towards the MC phase are smooth, i.e., there is a 
discontinuity in the second derivative of the current pro- 
file. On the other hand, the transition between LD and 
HD is an abrupt transition. The critical points a c and p c 
have the same dependency on the parameters k and 7. 
The transition line between LD and HD is a straight line 
in the a — P plane and is given by the condition a = p. 
Other works taking into account the biochemical cycle of 
ribosomes 27J find a more complicated relation between 
a c and p c , apparently due to another choice of the densi- 
ties in the reservoirs. However, here we obtain the same 
results with both the MCP and the mean-field analysis. 

For a fixed value of the translocation rate 7, the bound- 
aries between the MC and the LD and HD regimes are 
shifted with varying k, and the MC phase becomes larger 
than the MC region of standard ASEPs. Importantly, in 
the limiting case k — > 00 the critical points approach the 
values obtained for an ASEP. Since this limit corresponds 
to consider transitions rij = 1 — > tii = 2 occurring instan- 
taneously by neglecting the internal state of the particles, 
we therefore confirm that the results obtained with our 
model are consistent with previous findings. 

Figure [9] shows the outcomes of numerical simulations. 
As in the close-boundary case, the analytical results 
present deviations from the Monte Carlo simulations for 
low values of fc/7. Despite that, a study of the numerical 
phase diagram (Fig. |9j:) shows the same phenomenology 
of the analytical treatment of Fig. [8] Only the loca- 
tion of the critical points is inaccurate in the mean-field 
theory. For instance, notice that for some values of the 
parameters, in Fig. [9]d the mean-field predicts a LD-HD 



transition instead of the smooth LD-MC transition nu- 
merically found, i.e., the MC phase is reached for lower 
values of a and P (green circles). In other words, the 
numerically observed MC region is even larger than the 
one predicted by the mean-field approximation. 




(b) 



u 0.2 0.4 0.6 0.8 

a 

FIG. 9: (Color online). Monte Carlo simulations. Panels (a) 
and (b) show the current J as a function of a in (a) systems 
passing from LD to MC {fi = 1) and (b) from LD to HD 
phase (ft = 0.12). Full lines represent the theoretical predic- 
tions and circles are simulation points (k = 1 in black and 
k — 0.1 in green - light gray). The value of 7 is fixed to 1. 
Panel (c) shows the numerical phase diagram (different colors 
represent different densities p) for a system with k = 7 = 1. 
Drawing the phase diagram for other values of the parameters, 
we reproduce the features illustrated in Fig. [8] 



V. DISCUSSION AND CONCLUSIONS 

In this work we have proposed a model for mRNA 
translation based on an exclusion process. We have ex- 



tended previous models based on the ASEP (Asymmet- 
ric Simple Exclusion Process) by including the internal 
stepping cycle of the ribosomes, which corresponds to 
allowing the particles to have multiple internal states. 
The same model has been previously introduced in [35] 
to study the traffic of molecular motors. We have con- 
densed the whole biochemical cycle of the ribosome into 
two main steps: (i) finding the correct tRNA, which oc- 
curs with rate k, and (ii) translocation rate of the ribo- 
some to the next codon, which happens with rate 7. 

This extension is crucial in describing the underlying 
biological process, since previous ASEP-based models ne- 
glect that ribosomes can find and keep a correct tRNA 
during the waiting time due to the occupation of the next 
codon. 

The main result of this work is that the transitions 
among the different dynamical regimes of the system oc- 
cur at different critical points than the ones predicted by 
former ASEP models. These critical points depend on 
both fc and 7. For example, when fc/7 is small, the MC 
phase is substantially enlarged compared to a standard 
ASEP with hopping rate 7 and which ignores the inter- 
nal degree of freedom (A; — > 00) [42]. Crucially, this is 
the biologically relevant regime, as shown by estimates 
of the parameters based on experimental data. Thus, 
this model describes the biological system much more 
accurately, and its predictions can be readily validated 
by experimental measurements. 

The analysis of the system with periodic boundary con- 
ditions introduces the model and the formalism. This 
situation has been studied in |35j where the authors pro- 
pose a mean-field approach to analyse the case with peri- 
odic boundary conditions. Here we present a more com- 
plete approach which recovers the previous results and, 
in addition, makes possible the treatment of the model 
in the open boundary case. When the ratio between the 
transition rates k and 7 is high, the current profile be- 
comes symmetric and the value of p* moves toward the 
expected value of a standard exclusion process. The devi- 
ation between the mean-field and our simulations for low 
values of the transition rate k has yet to be understood. 
Furthermore, our model with periodic boundary condi- 
tions allows us to study whether there is a dominance of 
particles in state 1 or 2. In the biological system that 
we describe (even if in a coarse-grained perspective), a 
lattice with the majority of sites having n, = 2 repre- 
sents an mRNA in which ribosomes are carrying the cog- 
nate tRNAs and are waiting for hopping. This might 
be unfavourable when a finite number of charged tRNAs 
is available. Roughly speaking, in these conditions the 
charged tRNAs are kept by the ribosomes and cannot be 
used to translate other codons. In this sense the usage of 
resources is not optimized if a > A. This result suggests 
that under starvation or stress conditions, there might 
be a transition from the a > A to the a < A regime. 

We have shown that this extension of the model has im- 
portant consequences for the different boundary-induced 
transitions. Namely, depending on the ratio of k and 7, 



the sizes of the low density (LD), high density (HD) and 
maximal current (MC) phases in the a — /3 parameter 
space can change substantially, where a and /3 repre- 
sent respectively the initiation and termination rate of 
ribosomes. Crucially, the phase diagram coincides with 
the one obtained with the ASEP if fc/7 1, whereas 
if fc/7 — > 0, the maximal current phase is enlarged to a 
great extent, and the transitions from the LD and HD 
to the MC phase occur at much lower values of a and f3 
(depending on the value of fc, simulations show that the 
critical points are overestimated by the mean-field ap- 
proach and the MC region is even larger than predicted) . 

Based on experimental data, the translocation rate 7 
is estimated to be 7 = 35 s" 1 [37] (which is naturally as- 
sumed to be constant for each codon), and the ratio fc/7 
turns out to be in the range 0.05 — 3.38, depending on 
the codon (these values are estimates based on [38 ). For 
most codons (~ 87%) the ratio of the rates is smaller than 
1. Therefore, for physiological conditions, our model pre- 
dicts a much larger MC phase in the parameter space 
than previous ASEP models which completely neglect 
the internal state of the ribosomes or, equivalently, as- 
sume that fc/7 ^ I contrary to the biological conditions. 
In that unrealistic situation, the ribosomes would find 
the cognate tRNAs as soon as they translocate to the 
next codon. 

Thus one could naively think that the translation pro- 
cess is optimized to produce the largest possible number 
of proteins per unit time. In other words, the current J 
is maximized, and hence the system is in the MC phase. 
Although this assertion is not justified and there might be 
cases in which other effects prevent the translation rate 
to become maximal, such as particular configurations of 
slow codons downstream of the 5' end of the mRNA, 
competition for common resources and regulation at the 
level of translation, in this work we find signatures in 
this direction. In particular, the model predicts that the 
MC phase occupies a very large region in the parameter 
space. Furthermore, the density profiles experimentally 
observed [39] recall the MC density profile of Fig. [7)3. Fi- 
nally, ribosome recycling |43] might drive the system to 
lie in the region with the highest current. In fact, the 
injection parameter a can be decomposed into two com- 
ponents: one constant coefficient a being the affinity of 
freely diffusing ribosomes to bind to an open mRNA, and 
an increasing function of J which represents the proba- 
bility per unit time that a ribosome leaving the end of 
the mRNA will be recycled. It is clear that a might 
grow until the current balances the maximal current J* 
and as a result, a is larger or equal to the critical value 
a c . The ribosome recycling and its potential impact on 
translation regulation has been thoroughly investigated 
in @D]. 

Further studies on the model proposed in this work will 
address lattices with inhomogeneities (slow codons) and 
the influence of the size of the particles (ribosomes are 
known to cover around 9 codons). On the biological side, 
we are planning to implement simulations of real mRNA 
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sequences from the S. cerevisiae genome and validate the 
model with experiments at different levels. 
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